SAMHSA NSSATS Exploratory Data Analysis

years <- seq(2014, 2019) 
#years <- seq(2014, 2020)
forprofit_counts <- profit_states_1419 %>% 
  dplyr::select(contains("count_Forprofit_")) %>% 
  summarize(across(everything(), sum, na.rm = TRUE))

Warning: There was 1 warning in `summarize()`.
ℹ In argument: `across(everything(), sum, na.rm = TRUE)`.
Caused by warning:
! The `...` argument of `across()` is deprecated as of dplyr 1.1.0.
Supply arguments directly to `.fns` through an anonymous function instead.

  # Previously
  across(a:b, mean, na.rm = TRUE)

  # Now
  across(a:b, \(x) mean(x, na.rm = TRUE))
This warning is displayed once every 8 hours.
Call `lifecycle::last_lifecycle_warnings()` to see where this warning was generated.

nonprofit_counts <- profit_states_1419 %>% 
  dplyr::select(contains("count_Nonprofit_")) %>% 
  summarize(across(everything(), sum, na.rm = TRUE))

govt_counts <- profit_states_1419 %>% 
  dplyr::select(contains("count_Government_")) %>% 
  summarize(across(everything(), sum, na.rm = TRUE))


plot(years, forprofit_counts, type = 'l', main= "Treatment Facility Counts over time", ylim = c(1000, 8000), col = 'blue', xlab = 'Year', ylab = '# Treatment Centers')
lines(years, nonprofit_counts, type = 'l', main= "Non Profit Counts", col = 'red')

lines(years, govt_counts, type = 'l', main= "Government Counts", col = 'green')


legend("bottomleft", c("for profit", "non profit", "government"), col = c("blue", "red", "green"), pch ='-', lwd = 2)

library(ggplot2)
state_final_df %>%
  distinct(State_Name, year, .keep_all = TRUE) %>%
  group_by(year) %>%
  summarize(mean_drugusecalc = mean(druguse_calc),
            mean_aud = mean(aud_yr),
            mean_marijuana = mean(marijuana_month),
            mean_cocaine = mean(cocaine_yr),
            mean_heroin = mean(heroin_yr)) %>%
  mutate(prev_year= as.integer(year) - 1,
          year_range = paste(prev_year, "-", year)) %>%
   ungroup() %>%
   filter(year <= 2018) %>% #interested in only data until 2018
   ggplot(., aes(x = year_range, group = 1)) + 
   geom_line(aes(y = mean_cocaine, color = "Cocaine")) + 
   geom_line(aes(y = mean_heroin, color = "Heroin")) + 
   geom_line(aes(y = mean_aud, color = "Alcohol Use Disorder")) + 
   geom_line(aes(y=mean_marijuana, color = "Marijuana")) + 
   scale_color_manual(values = c("Cocaine" = "red", "Heroin" = "blue", "Alcohol Use Disorder" = "purple", "Marijuana" = "green")) +
   labs(y = "Average Prevalence in 50 States + DC", title = "NSDUH Alcohol Use Disorder and Drug Use Prevalence over Time")# + 

   #scale_x_discrete(labels = year_range)
  # #labs(y = "Average Prevalence in 50 States + DC", title = "NSDUH Alcohol Use Disorder and Marijuana Use Prevalence over Time")
  # #labs(y = "Average Prevalence in 50 States + DC", title = "NSDUH Cocaine and Heroin Use Prevalence over Time")

state_final_df %>%
  #distinct(State_Name, year, .keep_all = TRUE) %>%
  filter(year >= 2015) %>%
  group_by(year, profit_type) %>%
  summarize(num_centers = sum(count_tx, na.rm = TRUE)) %>%
  ungroup() %>%
  pivot_wider(id_cols = year, names_from = profit_type, values_from = num_centers) %>%
  ggplot(., aes(x = year, group = 1)) + 
  geom_line(aes(y = Forprofit, color = "For-Profit Centers")) + 
  geom_line(aes(y = Nonprofit, color = "Non-Profit Centers")) + 
  geom_line(aes(y = Government, color = "Government Centers")) + 
  scale_color_manual(values = c("For-Profit Centers" = "green", "Non-Profit Centers" = "blue", "Government Centers" = "red")) +
  labs(y = "Sum Treatment Centers in 50 States + DC", title = "Treatment Center Counts by Profit Status over Time")

`summarise()` has grouped output by 'year'. You can override using the `.groups` argument.

using map tutorial from: https://jtr13.github.io/cc19/different-ways-of-plotting-u-s-map-in-r.html#using-ggplot2-package state abbreviations from: https://www.ssa.gov/international/coc-docs/states.html

library(ggplot2)
library(maps)
library(mapdata)
library(plotly)

Registered S3 method overwritten by 'data.table':
  method           from
  print.data.table     
Registered S3 method overwritten by 'htmlwidgets':
  method           from         
  print.htmlwidget tools:rstudio

Attaching package: ‘plotly’

The following object is masked from ‘package:ggplot2’:

    last_plot

The following object is masked from ‘package:stats’:

    filter

The following object is masked from ‘package:graphics’:

    layout

library(sf)

Linking to GEOS 3.11.0, GDAL 3.5.3, PROJ 9.1.0; sf_use_s2() is TRUE

usa <- map_data('usa')
state <- map_data("state")

#merge state geo file to the for-profit information
state <- state %>% 
              left_join(state_abbrevs_df, by = c("region" = "state_names")) %>% 
              left_join(profit_states_1419, by = c('state_abbrevs' = 'STATE')) 
                                      

#plot the number of for-profit treatment centers from 2015-2019
for (i in seq(2015, 2019))
{
  plot_forprofits <- ggplot(data=state, aes(x=long, y=lat, fill= !!sym(paste0("count_Forprofit_", i, sep = "")), group=group)) + 
                      scale_fill_gradient(low = "grey", high = "green", limits = c(0,800)) +
                      geom_polygon(color = "white") + 
                      ggtitle(paste('# For-Profit Treatment Centers in', i)) + 
                      coord_fixed(1.3)
  print(plot_forprofits)
}

ggplot(data=state, aes(x=long, y=lat, fill= change_Forprofit_1519, group=group)) + 
                      scale_fill_gradient(low = "grey", high = "green", limits = c(0,400)) +
                      geom_polygon(color = "white") + 
                      ggtitle('Change in # For-Profit Treatment Centers from 2015 - 2019') + 
                      coord_fixed(1.3)

for (i in seq(2015, 2019))
{
  plot_nonprofits <- ggplot(data=state, aes(x=long, y=lat, fill= !!sym(paste0("count_Nonprofit_", i, sep = "")), group=group)) + 
                      scale_fill_gradient(low = "grey", high = "blue", limits = c(0,900)) +
                      geom_polygon(color = "white") + 
                      ggtitle(paste('# Non-Profit Treatment Centers in', i)) + 
                      coord_fixed(1.3)
  print(plot_nonprofits)
}

ggplot(data=state, aes(x=long, y=lat, fill= change_Nonprofit_1519, group=group)) + 
                      scale_fill_gradient(low = "grey", high = "blue", limits = c(0,100)) +
                      geom_polygon(color = "white") + 
                      ggtitle('Change in # Non-Profit Treatment Centers from 2015 - 2019') + 
                      coord_fixed(1.3)

for (i in seq(2015, 2019))
{
  plot_govt <- ggplot(data=state, aes(x=long, y=lat, fill= !!sym(paste0("count_Government_", i, sep = "")), group=group)) + 
                      scale_fill_gradient(low = "grey", high = "red", limits = c(0,250)) +
                      geom_polygon(color = "white") + 
                      ggtitle(paste('# Government Treatment Centers in', i)) + 
                      coord_fixed(1.3)
  print(plot_govt)
}

ggplot(data=state, aes(x=long, y=lat, fill= change_Government_1519, group=group)) + 
                      scale_fill_gradient(low = "grey", high = "red", limits = c(0,50)) +
                      geom_polygon(color = "white") + 
                      ggtitle('Change in # Government Treatment Centers from 2015 - 2019') + 
                      coord_fixed(1.3)

# state_name <- "California"
# 
# par(mar=c(5, 4, 4, 12), xpd=TRUE)
# 
# plot(yearly_nsduh_data[yearly_nsduh_data$stname == state & yearly_nsduh_data$outcome == "UDPYILAL", ]$start_year, yearly_nsduh_data[yearly_nsduh_data$stname == state & yearly_nsduh_data$outcome == "UDPYILAL", ]$est_total, type = 'o', ylim= c(0,10000), main = paste(state, "Substance Use Measures Over Time") ,xlab = "year", ylab = "Thousands of People in 2015-18")
# 
# lines(yearly_nsduh_data[yearly_nsduh_data$stname == state & yearly_nsduh_data$outcome == "ILLEMMON", ]$start_year, yearly_nsduh_data[yearly_nsduh_data$stname == state & yearly_nsduh_data$outcome == "ILLEMMON", ]$est_total, col = 'blue', type = 'o')
# 
# lines(yearly_nsduh_data[yearly_nsduh_data$stname == state & yearly_nsduh_data$outcome == "BNGDRK", ]$start_year, yearly_nsduh_data[yearly_nsduh_data$stname == state & yearly_nsduh_data$outcome == "BNGDRK", ]$est_total, col = 'red', type = 'o')
# 
# lines(yearly_nsduh_data[yearly_nsduh_data$stname == state & yearly_nsduh_data$outcome == "TXNPILAL", ]$start_year, yearly_nsduh_data[yearly_nsduh_data$stname == state & yearly_nsduh_data$outcome == "TXNPILAL", ]$est_total, col = 'green', type = 'o')
# 
# lines(yearly_nsduh_data[yearly_nsduh_data$stname == state & yearly_nsduh_data$outcome == "PNRNMYR", ]$start_year, yearly_nsduh_data[yearly_nsduh_data$stname == state & yearly_nsduh_data$outcome == "PNRNMYR", ]$est_total, col = 'purple', type = 'o')
# 
# legend('topright',legend = c("Substance Use Disorder", "Drug use in past month", "Binge alcohol use", "Needing, not receiving tx", "pain reliever misuse"), col = c("black", "blue", "red", "green", "purple"), lty = 1,  inset=c(-0.5, 0))

plot(state_final_df[state_final_df$profit_type == "Forprofit", ]$sud, state_final_df[state_final_df$profit_type == "Forprofit",]$tx_percapita, xlab = "Substance Use Disorder Counts (in thousands)", ylab = "For profit treatment centers per capita")
abline(lm(tx_percapita ~ sud, data = state_final_df[state_final_df$profit_type == "Forprofit", ]))

plot(state_final_df[state_final_df$profit_type == "Nonprofit", ]$sud, state_final_df[state_final_df$profit_type == "Nonprofit",]$tx_percapita, col = "blue", xlab = "Substance Use Disorder counts (in thousands) ", ylab = "Nonprofit treatment centers per capita")
abline(lm(tx_percapita ~ sud, data = state_final_df[state_final_df$profit_type == "Nonprofit", ]))

plot(state_final_df[state_final_df$profit_type == "Government", ]$sud, state_final_df[state_final_df$profit_type == "Government",]$tx_percapita, col = "green", xlab = "Substance Use Disorder counts (in thousands)", ylab = "Government treatment centers per capita")
abline(lm(tx_percapita ~ sud, data = state_final_df[state_final_df$profit_type == "Government", ]))

NA
NA

Look at the distribution of each input variable of interest at each point in time

create_hist <- function(df, year_param, profit_param, col_param)
{
  subset_data <- df %>%
  filter(year == year_param, profit_type == profit_param)
  hist(subset_data[[col_param]], main = paste("Histogram of ",year_param, profit_param, col_param), xlab = col_param)
}

#substance use measures are roughly normal
create_hist(state_final_df, 2016, "Forprofit", "binge_drinking")

create_hist(state_final_df, 2017, "Forprofit", "binge_drinking")

create_hist(state_final_df, 2018, "Forprofit", "binge_drinking")

create_hist(state_final_df, 2019, "Forprofit", "binge_drinking")

create_hist(state_final_df, 2016, "Forprofit", "illicit_druguse")

create_hist(state_final_df, 2017, "Forprofit", "illicit_druguse")

create_hist(state_final_df, 2018, "Forprofit", "illicit_druguse")

create_hist(state_final_df, 2019, "Forprofit", "illicit_druguse")

create_hist(state_final_df, 2016, "Forprofit", "needing_tx")

create_hist(state_final_df, 2017, "Forprofit", "needing_tx")

create_hist(state_final_df, 2018, "Forprofit", "needing_tx")

create_hist(state_final_df, 2019, "Forprofit", "needing_tx")

create_hist(state_final_df, 2016, "Forprofit", "sud")

create_hist(state_final_df, 2017, "Forprofit", "sud")

create_hist(state_final_df, 2018, "Forprofit", "sud")

create_hist(state_final_df, 2019, "Forprofit", "sud")

#for profit treatment centers are heavily skewed for all years
create_hist(state_final_df, 2016, "Forprofit", "count_tx")

create_hist(state_final_df, 2017, "Forprofit", "count_tx")

create_hist(state_final_df, 2018, "Forprofit", "count_tx")

create_hist(state_final_df, 2019, "Forprofit", "count_tx")

#non profit treatment centers are heavily skewed for all years
create_hist(state_final_df, 2016, "Nonprofit", "count_tx")

create_hist(state_final_df, 2017, "Nonprofit", "count_tx")

create_hist(state_final_df, 2018, "Nonprofit", "count_tx")

create_hist(state_final_df, 2019, "Nonprofit", "count_tx")

#government treatment centers are heavily skewed for all years
create_hist(state_final_df, 2016, "Government", "count_tx")

create_hist(state_final_df, 2017, "Government", "count_tx")

create_hist(state_final_df, 2018, "Government", "count_tx")

create_hist(state_final_df, 2019, "Government", "count_tx")

Investigating treatment centers per capita over time and by profit status

forprofit_txpercapita <- state_final_df %>% 
  filter(profit_type == "Forprofit") %>% 
  group_by(year) %>%
  summarize(avg_tx_percapita = mean(tx_percapita, na.rm = TRUE),
            min_tx_percapita = min(tx_percapita, na.rm = TRUE),
            max_tx_percapita = max(tx_percapita, na.rm = TRUE))
  
nonprofit_txpercapita <- state_final_df %>% 
  filter(profit_type == "Nonprofit") %>% 
  group_by(year) %>%
  summarize(avg_tx_percapita = mean(tx_percapita, na.rm = TRUE),
            min_tx_percapita = min(tx_percapita, na.rm = TRUE),
            max_tx_percapita = max(tx_percapita, na.rm = TRUE))

govt_txpercapita <- state_final_df %>% 
  filter(profit_type == "Government") %>% 
  group_by(year) %>%
  summarize(avg_tx_percapita = mean(tx_percapita, na.rm = TRUE),
            min_tx_percapita = min(tx_percapita, na.rm = TRUE),
            max_tx_percapita = max(tx_percapita, na.rm = TRUE))



plot(forprofit_txpercapita$year, forprofit_txpercapita$avg_tx_percapita, type = 'o', col = 'blue', ylim = c(0.000001, 0.00005))
#lines(forprofit_txpercapita$year, forprofit_txpercapita$min_tx_percapita, type = 'l', col = 'blue', lty = 2)
#lines(forprofit_txpercapita$year, forprofit_txpercapita$max_tx_percapita, type = 'l', col = 'blue', lty = 2)

lines(nonprofit_txpercapita$year, nonprofit_txpercapita$avg_tx_percapita, type = 'o', col = 'red' )

#lines(nonprofit_txpercapita$year, nonprofit_txpercapita$min_tx_percapita, type = 'l', col = 'red', lty = 2 )
#lines(nonprofit_txpercapita$year, nonprofit_txpercapita$max_tx_percapita, type = 'l', col = 'red', lty = 2 )

lines(govt_txpercapita$year, govt_txpercapita$avg_tx_percapita, type = 'o', col = 'green')

#lines(govt_txpercapita$year, govt_txpercapita$min_tx_percapita, type = 'l', col = 'green', lty = 2)
#lines(govt_txpercapita$year, govt_txpercapita$max_tx_percapita, type = 'l', col = 'green', lty = 2)

state_final_df[is.na(state_final_df$tx_percapita),]

#Connecticut, South Carolina, Ohio, Mississippi, New York have the fewest for-profit treatment centers per capita
#Maine, Idaho, Utah, North Dakota, Kentucky have the most for-profit treatment centers per capita 
state_final_df %>%
  group_by(State_Abbrev, profit_type) %>%
  summarize(avg_txpercapita = mean(tx_percapita, na.rm = TRUE)) %>%
  filter(profit_type == "Forprofit") %>%
  arrange(avg_txpercapita) %>%
  ungroup() %>%
  left_join(state, by = c("State_Abbrev" = "state_abbrevs")) %>%
ggplot(aes(x=long, y=lat, fill= avg_txpercapita, group=group)) + 
                      scale_fill_gradient(low = "grey", high = "green", limits = c(0,0.0001)) +
                      geom_polygon(color = "white") + 
                      ggtitle('Average Nonprofit Treatment Centers per capita 2015-2019') + 
                      coord_fixed(1.3)

`summarise()` has grouped output by 'State_Abbrev'. You can override using the `.groups` argument.

#Virginia, Texas, Idaho, South Carolina, Georgia have the fewest non-profit treatment centers per capita
#Hawaii, Alaska, Maine, Vermont, Wyoming have the most non-profit treatment centers per capita 
state_final_df %>%
  group_by(State_Abbrev, profit_type) %>%
  summarize(avg_txpercapita = mean(tx_percapita, na.rm = TRUE)) %>%
  filter(profit_type == "Nonprofit") %>%
  arrange(avg_txpercapita) %>%
  ungroup() %>%
  left_join(state, by = c("State_Abbrev" = "state_abbrevs")) %>%
ggplot(aes(x=long, y=lat, fill= avg_txpercapita, group=group)) + 
                      scale_fill_gradient(low = "grey", high = "blue", limits = c(0,0.0001)) +
                      geom_polygon(color = "white") + 
                      ggtitle('Average Nonprofit Treatment Centers per capita 2015-2019') + 
                      coord_fixed(1.3)

`summarise()` has grouped output by 'State_Abbrev'. You can override using the `.groups` argument.

#New Hampshire, Tennessee, Pennsylvania, Delaware, Kentucky have the fewest government treatment centers per capita
#Alaska, North Dakota, South Dakota, Wyoming, New Mexico have the most government treatment centers per capita 
state_final_df %>%
  group_by(State_Abbrev, profit_type) %>%
  summarize(avg_txpercapita = mean(tx_percapita, na.rm = TRUE)) %>%
  filter(profit_type == "Government") %>%
  arrange(avg_txpercapita) %>%
  ungroup() %>%
  left_join(state, by = c("State_Abbrev" = "state_abbrevs")) %>%
ggplot(aes(x=long, y=lat, fill= avg_txpercapita, group=group)) + 
                      scale_fill_gradient(low = "grey", high = "red", limits = c(0,0.00005)) +
                      geom_polygon(color = "white") + 
                      ggtitle('Average government Treatment Centers per capita 2015-2019') + 
                      coord_fixed(1.3)

`summarise()` has grouped output by 'State_Abbrev'. You can override using the `.groups` argument.

Create a csv with the averages per state between 2015 - 2019 and write to csv

state_final_df %>%
  filter(year >= 2015) %>% #model will only use 2015 - 2019
  group_by(State_Name, year) %>%
  pivot_wider(id_cols = c(State_Name, year), names_from = profit_type, values_from = count_tx) %>%
  ungroup() %>%
  arrange(year) %>%
  group_by(State_Name) %>%
  summarize(Forprofit_Count2015 = first(Forprofit),
            Forprofit_Change1519 = last(Forprofit) - first(Forprofit),
            Nonprofit_Count2015 = first(Nonprofit),
            Nonprofit_Change1519 = last(Nonprofit) - first(Nonprofit),
            Government_Count2015 = first(Government), 
            Government_Change1519 = last(Government) - first(Government)
) %>%
  arrange(desc(Forprofit_Change1519)) %>%
  write.csv("State_Profit_Counts.csv")
  #mutate(increase = ifelse())

Create a csv of the drug use by state

state_final_df %>%
  filter(year >= 2015) %>% #model will only use 2015 - 2019
  group_by(State_Name) %>%
  arrange(year) %>%
   summarize(Druguse2013_14 = first(prev2years_drugusecalc),
             Cocaine2013_14 = first(prev2years_cocaine),
             Heroin2013_14 = first(prev2years_heroin),
             AUD2013_14 = first(prev2years_aud),
             Marijuana2013_14 = first(prev2years_marijuana))%>%
#             Forprofit_Change1519 = last(Forprofit) - first(Forprofit),
#             Nonprofit_Count2015 = first(Nonprofit),
#             Nonprofit_Change1519 = last(Nonprofit) - first(Nonprofit),
#             Government_Count2015 = first(Government), 
#             Government_Change1519 = last(Government) - first(Government)
   arrange(desc(Druguse2013_14)) %>%
    write.csv("State_Druguse_Counts.csv")

state_yearly_averages<- state_final_df %>%
  distinct(State_Abbrev, profit_type, year, .keep_all = TRUE) %>%
  group_by(year) %>%
  summarize(avg_druguse = mean(illicit_druguse),
         avg_drinking = mean(binge_drinking),
         avg_od = mean(RATE)) %>%
  arrange(year)

plot(state_yearly_averages$year, state_yearly_averages$avg_druguse, type = 'l')

plot(state_yearly_averages$year, state_yearly_averages$avg_drinking, type = 'l')

plot(state_yearly_averages$avg_od, state_yearly_averages$avg_od, type = 'l')

state_final_df %>%
  select(State_Name, profit_type, year, count_tx, cocaine_yr, marijuana_month, aud_yr, heroin_yr, druguse_calc) %>%
  #group_by(State_Name, profit_type) #%>%
  pivot_wider(id_cols = c(State_Name, profit_type), names_from = year, values_from = count_tx, unused_fn = first) #%>%
  #mutate(average_1419 = mean(select("2014", "2015", "2016", "2017", "2018", "2019")))#%>%
  #pivot_wider(id_cols = c(State_Name, year), names_from = profit_type, values_from = count_tx, unused_fn = first) #%>%
  #pivot_wider(id_cols = c(State_Name), names_from = profit_type, values_from=c("2014", "2015", "2016", "2017", "2018", "2019"), unused_fn = first)

print("Test")